rm(list=ls())
library(sf)
library(ggplot2)
library(RColorBrewer)
library(ggthemes)
library(tidycensus)
library(tidyverse)
library(viridis)

key_census <- ""
dir <- ""
#read in boundaries 
hs_boundaries <- st_read( paste(dir, "zone_boundaries_zoc_included.shp", sep="") ) 
# subset to ZOC 
zoneids <- c(26, 62, 74, 50, 76, 77, 12, 73, 
             55, 84, 82, 38, 81, 82, 84, 
             83, 58, 20, 37 )
zone_boundaries <- hs_boundaries
zone_boundaries["zone"] <- 0
zone_boundaries$zone[which(zone_boundaries$OBJECTID %in% zoneids)] <- 1
zone_boundaries <- subset(zone_boundaries, zone_boundaries$zone==1)
names(zone_boundaries$OBJECTID) <- "zoneid"  
colnames(zone_boundaries)[colnames(zone_boundaries)=="OBJECTID"] <- "zoneid"

# read in ACS data 
tarr <- get_acs(geography = "tract", year=2010, variables = "B19013_001",state = "CA", county = "Los Angeles", geometry = TRUE)
tarr<- st_transform(tarr, 4326)

# merge blockids with zone information 
blocks_with_hs <- st_join(tarr, join=st_intersects, left=TRUE, hs_boundaries["OBJECTID"])
blocks_with_zones <- st_join(blocks_with_hs, join=st_intersects, left=TRUE, zone_boundaries["zoneid"])
hs <- subset(blocks_with_hs, !is.na(blocks_with_hs$OBJECTID))
income_cuts <- c(0, 35223.25, 47220, 70636, 232935)
hs["income_quartiles"] <- cut(hs$estimate, breaks=income_cuts, include.lowest=FALSE)
hs$quartile <- NA
hs[which(hs$estimate <=35223.25), ]$quartile <- 1
hs[which(hs$estimate >35223.25 & hs$estimate<=47220 ), ]$quartile <- 2
hs[which(hs$estimate >47220 & hs$estimate<=70636 ), ]$quartile <- 3
hs[which(hs$estimate >70636 ), ]$quartile <- 4
#scale_fill_gradient(low = "yellow", high = "red", aesthetics = "fill") 
  

blues_palette <- rev(brewer.pal(n = 4, name = "Blues"))
ggplot() +
  geom_sf(data = hs, aes(fill = as.factor(quartile)), size=0.3) +
  geom_sf(data = hs_boundaries, fill = NA, colour = "black", size = .6) +
  geom_sf(data = zone_boundaries, fill = NA, colour = "red", size = .6) +
  theme_bw() +
  scale_fill_manual(
    "Tract median income quartile (2010)",
    values = blues_palette,
    labels = c("4th Quartile", "3rd Quartile", "2nd Quartile", "1st Quartile")
  ) +
  theme_map()

ggsave(paste(dir,"censustractincome.pdf", sep=""))
